function [syncM]=SLopt(dataM,fsample,filterLF,filterHF)
    
    %SL Parameter
        

    epochLength=size(dataM,2);
    
	[param.lag,param.embDim,param.theiler,param.pRef,param.speed,param.window]=autoSetParameters(fsample,filterLF,filterHF,epochLength);
%     param.lag = 10;
%     param.embDim = 10;
%     param.theiler = 100;
% 	param.pRef = 0.01;
%     param.speed = 16;
% 	param.window = 410;
    param
    [syncM,avgSync,avgSyncXCh]= multiSL(param,dataM);

    
    
    function [lag,embDim,theiler,pRef,speed,window]=autoSetParameters(fsample,filterLF,filterHF,epochLength)
    
    lag=round(fsample/(4*filterHF));
    if(lag<1)
        lag=1;
    end
    embDim=1+round(fsample/(filterLF*lag));
    if(embDim>20)
        embDim=20;
    end
    theiler=lag*embDim;
    pRef=0.01;
    speed=16;
    window=round(0.1*epochLength);
    
    while((theiler*2)>window)
        lag=lag-1;
        theiler=lag*embDim;
    end
    
    if(lag<1)
        lag=1;
    end

    
    
function [distance]=computeVectorDistance(i, j, k,param,dataM)
    distance = 0;
 
    for  a=0:param.embDim-1
        %distance += Math.pow((myData.getData()[k][i+a*lag]-myData.getData()[k][j+a*lag]),2);
        distance= distance + (dataM(k,i+a*param.lag)-dataM(k,j+a*param.lag))^2;
    end
    distance = sqrt(distance);

    
    
     
	
function [cutoff]=computeCutoff(i,k,param,dataM)
	
		%double[] distances = new double[4097];
    nSamples=size(dataM,2);
    epochLength=nSamples;
    distances=zeros(1,nSamples);%%4097 check if important??
    tel = 0;
    hulp = 0; 
    cutoff = 0;
		
		%for (int  a=1; a<myData.getEpochLength()-lag*embDim; a++)
    for a=1:(epochLength-param.lag*param.embDim)-1
        if ((abs(a-i)>param.theiler) && (abs(a-i)<param.window))			
            tel=tel+1;
            distances(tel) = computeVectorDistance(i,a,k,param,dataM);
        end
    end
    
    d2=sort(distances(1:tel));
    indx=round(param.pRef*tel);
    if(indx~=0)
        cutoff=d2(indx);
    else
        cutoff=0;
    end
%     for a=1:tel-1
%         for b=1:tel-1
%             if (distances(b+1) < distances(b))
%                 hulp = distances(b+1);
%                 distances(b+1) = distances(b);
%                 distances(b) = hulp;
%             end
%         end
%     end
%     
%     cutoff = distances(round(param.pRef*tel));
		
	
	
function [syncM,avgSync,avgSyncXCh]= multiSL(param,dataM)
	
    nCh=size(dataM,1);
    nSamples=size(dataM,2);
    syncM=zeros(nCh,nCh);
    avgSync=0;
    avgSyncXCh=zeros(1,nCh);
    
    lastChannel=nCh;
    epochLength=nSamples;

    rCutoff = zeros(1,nCh);
    chanHits = zeros(1,nSamples);
    hit = zeros(1,nCh);

    length = epochLength - param.lag*param.embDim;

    %for (int i=1; i<length; i+=speed)
    for i=1:param.speed:length-1
        %compute cutoff distance for time i, for each channel k
        for k=1:lastChannel
            rCutoff(k) = computeCutoff(i,k,param,dataM);
        end

        for j=1:length-1
            if ((abs(i-j)>param.theiler) && (abs(i-j)<param.window))	

                for k=1:lastChannel					
                    hit(k) = 0;
                    if (computeVectorDistance(i,j,k,param,dataM) < rCutoff(k))
                        hit(k) = 1;
                        chanHits(k)=chanHits(k)+1;
                    end %if computeVectorDistance
                end %for k

                for k1=1:lastChannel
                    for k2=1:lastChannel
                        if (hit(k1) && hit(k2) && (k1~=k2)) 			
                            syncM(k1,k2)=syncM(k1,k2)+ 1;
                        end % if hit hit            
                    end  %for k2
                end %for k1

            end %if theiler window
        end%for length
    end%for speed
		
		
    totCount = 0;
    for  k1=1:lastChannel
        chanCount = 0;
        for k2=1:lastChannel
            if ( k1~=k2)
                chanCount = chanCount+1;
                totCount = totCount+1;
                if ((chanHits(k1)+chanHits(k2)) > 0) 
                    syncM(k1,k2) = syncM(k1,k2)/ (0.5*chanHits(k1)+0.5*chanHits(k2));
                    avgSync = avgSync +syncM(k1,k2);
                    avgSyncXCh(k1) = avgSyncXCh(k1) + syncM(k1,k2);
                end %if chanHits
            
            end%if k1!=k2
        
        end% for k2
      
        if (chanCount>0) 
            avgSyncXCh(k1) = avgSyncXCh(k1)/chanCount;
        end
        
    end%for k1
    
    if (totCount>0) 
        avgSync = avgSync/totCount;
    end
		

	

